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Abstract 

In this article a stochastic cellular automata model is examined, which 
has been developed to study a "small" world, where local changes may 
noticeably alter global characteristics. This is applied to a climate model, 
where global temperature is determined by an interplay between atmo- 
spheric carbon dioxide and carbon stored by plant life. The latter can be 
relased by forest fires, giving rise to significant changes of global conditions 
within short time. 

Keywords: Cellular Automata; Forest Fire; Climate Dynamics; Carbon 
Cycle 

PACS Nos.: 64. 60. ah, 89.75.-k, 89.75.Fb, 89.75.kd 

1 Introduction 

1.1 A Brief Note on Cellular Automata 

Cellular Automata (CA), as studied first by von Neumann and Ulam pQ, provide 
systems both interesting from a fundamental point of view and for practical 
purposes of model builders. Among the systems most frequently studied are 
CA versions of daisyworld [HOE] and forest fire models EE M ■ 

The main idea is to study a grid or network of identical units (cells), which 
each can take discrete states, and where the time evolution of each cell is gov- 
erned by simple local rules: The state Si(t) of a cell labelled by i at timestep t 



is only determined by its previous state Sj(t — 1) and the states Sj k (t — 1) of a 
small number k of neighbours, labelled by jk- 

There arc different ways how to precisely define a neighbourhood, but typ- 
ically the qualitative behaviour is typically unaffected by the details [TO] [IT] . 
This can be seen as a footprint of universality - the fact that microscopic de- 
tails do not influence features of systems which are "critical" in the sense of 
Wilson [T^j . The fact that all rules are local implies a finite maximum speed for 
the propagation of all effects. 

Originally the rules governing cellular automata were strictly deterministic. 
One can extend this formalism (and make it both more flexible and more real- 
istic) when including stochastic elements, i.e. give probabilities how the state 
of a cell depends on its own previous state and those of its neighbours. On can 
also replace the few discrete states by a (quasi)continuum to gain additional 
flexibility. The system studied in this article belongs to this class of Stochas- 
tic Cellular Automata (SCA) models, which the authors apply to a particular 
feature of climate-relevant dynamics. 

1.2 A Brief Note on the Carbon Cycle 

It is well-known that the global temperature depends strongly on the amount 
of greenhouse gases contained in the atmosphere. While H2O vapor is the 
most important of them, its dynamic is extremely complicated and hard to 
treat reliably in models. Since the amount of H 2 vapor strongly depends on 
temperature, the water vapor cycle may largely act as a positive feedback cycle, 
enhancing the effects of other greenhouse gases. 

For these reasons the authors focus on CO2 , being the most important other 
greenhouse gas. Atmospheric CO2 is crucial for plant life and one typically 
expects enhanced growth for larger amount of CO2 available. This is indeed 
sometimes used as an argument against global warming [13] . 

But in particular when accompanied by climate changes, enhanced plant 
growth can also mean increased vulnerability to large-scale forest fire, which can 
release considerable amounts of CO2 in very short time [14] (though balancing 
aspects seem possible [15]). 

There exists a number of realistic climate models, which attempt to directly 
describe the terrestric system [16] . These models, which are by now quite suc- 
cessful, are mostly based on the (simplified) treatment of atmosphere and ocean 
dynamics, but mostly neglect the plant life Carbon cyle. Therefore the present 
model, though significantly more simplistic, can be regarded as complementary 
to these approaches. 
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2 The Model 



2.1 The Basic Idea 

In this model one studies patches of land, which are, in the course of time, 
inhabited by plant life. During their growth process, these plants can store 
carbon from the atmosphere. The growth rate depends both on temperature 
(determined by geographical latitude and carbon dioxide in the atmosphere) and 
the amount of carbon available (i.e. again carbon dioxide in the atmosphere). 

For each cell, optimum growth temperature and temperature tolerance are 
first set randomly and then subject to a simple evolutionary process. Uninhab- 
ited cells can become inhabited by spreading from neighbours. Their character- 
istic parameters are determined as a weighted average of those of all neighbours, 
affected by mutations. 

In addition, there arc spontaneous outbreaks of forest fires, and patches are 
more vulnerable to them, if they are more heavily overgrown. The fires can 
spread, and burning cells immediately release all carbon stored in them to the 
atmosphere in form of carbon dioxide. 

Remark: Since the amount of oxygen in the atmosphere is not noticably 
changed during all processes, in the following the distinction between "carbon" 
and "carbon dioxide" will be dismissed and always the word "carbon" is used, 
where it is implicitly assumed that a sufficient amount of oxygen is available 
for all processes and that carbon is deposited in the atmosphere in the form of 
carbon dioxide. 

In the process of biological degradation and during fires, not only CO2, but 
also considerable amounts of H 2 are released. While, for reasons discussed in 
section [L2| we do not attempt to explicitly include water vapor effects, they may 
be - at least at some timescales - partially covered too. 

2.2 Relation to Previous Work 

The system studied in this approch is closely related to the Clar-Drossel-Schwabl 
forest-fire model [3 [8] which also combines slow growth with rapid burning of 
clusters. In contrast to the Clar-Drossel-Schwabl model, the present approach 
features stochastic propagation of the fire and - more important - the coupling 
to a global variable, to be interpreted as the atmospheric carbon dioxide content 
of the model world. 

2.3 Implementation 

The model is set up on a rectangular grid, described by several [iV, x Nj]- 
matrices, which are summarized and explained in table[l] Each cell represents a 
patch of land, possibly inhabited by a certain amount of plant life with individual 
characteristics. Boundary conditions are chosen periodic in the first (i) and open 
in the second (j) direction. 
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matrix 


range 


meaning 


L 


% G {0, 1} 


life/inhabitance Matrix 


C 


Cij G [0, 1] 


amount of carbon stored in a specific cell 


T 


Tij e R 


optimal temperature for a cell 


S 




temperature tolerance (width of distribution) 


B 


&« G {0, 1} 


flag for burning cell (used only in forest fire phase) 



Figure 1: System matrices 



A simulation run is performed as follows: 

Initialization: Cells are inhabited with starting density p s tart, the parame- 
ters of inhabited cells are determined as = U, Tij — U and <zy = f +fi-(2U— 1) 
with U taken individually from a uniform [0, l)-distribution. Now the following 
cycle is repeated iV s teps times: 

1. Spreading: Each cell with at least one inhabited neighbour has a chance 
to become inhabited as well. The probability is given bjFI 

i'j' 

where the prime at the sum indicates summation over next neighbours, 
with diagonal elements (i ± 1, j ± 1) weighted with the Moore factor [TU] 
ttm = \ ■ In case the cell becomes inhabited, its parameters are determined 
as 

T y = ^ E' c ^'^'+A*-(2E/-l), (2) 

i'j' 

u » = Jf £' c ^' T *v 0- + m • - I)) , (3) 

i'j' 

with N :— ^2 Ci'j' and U individually taken from a uniform random 
distribution on [0, 1). The row-averaged values 

lijTij and <jj = =-^- yltijo-ij (4) 



3 p / j 3 *3 t *"- u "3 ~ sr^ y 'J 

are stored for later time series analysis. 



In principle, one could introduce an additional multiplicative parameter « apre ad S (0, 1] in 
front of the sum, the same way it is done for the following processes. Our choice of « sprea( j = 1 
effectively fixes the timescale. Accordingly the rates of all other processes (except spreading 
of forest fires) are determined relative to this spreading rate. 
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2. Determining global variables (1): The carbon content C of the atmo- 
sphere is determined as 



C= 1- 



NiNj 

and the temperature is determined as 

j - 1 

2} = Q!Tvar — 1 + (1 - ttTvar) C . (6) 

The value of C is stored for later time-series analysis. 
3. Growth: Each cell growths with its individual rate, 

dj(t + 1) = medjo, Cij(t) + ^growth , Tij, 

Pij), l} (7) 

(where "med" denotes the median value) with the temperature adaptation 
model function 



2a,, 



(p(Tj, t^, (Tij) = Nexp \ \ - fthosta > (8) 

7 ij 



where the normalization N = Nfc;, oy) is chosen such that 

V?(Tj, Ty, (Tij) = 1 - /ihostil ■ (9) 



This simple Gaussian shape already shows the most important features 
one would expect of such a function. It allows both extreme specialization 
to a certain temperature and adaptation to a wider range. For /ihostil > 
there may be temperature regions where growth is negative. If this causes 
Cij(t + 1) to be = 0, the cell dies. 

4. Determining global variables (2): Again Carbon content and temper- 
ature are determined according to equations ([5| and 

5. Forest fires: Each inhabited cell (i, j) has a chance to catch fire, with 
probability 

p = a &Te T j c i j. (10) 

If there any burning cells (cells that caught fire) , a complementary simu- 
lation on a short timescale is started. The fire may spread from burning 
cells (bij = 1) to their neighbours: 

(a) In each step of the burning process, each inhabited cell can catch fire 
with probability 

i'ji 
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param. 


range 


meaning 


Ni 


eN> 2 


linear size of the grid in first direction (open b.c.) 


Nj 


eN> 2 


linear size of the grid in second direction (periodic b.c.) 


^max 


e2 N 


simulation time 


Pstart 


G (0, 1] 


initial density of inhabited cells 


afire 


G M+ 


probability of a cell spontaneously catching fire 


^growth 


G M+ 


measure for growth rate 


^■hostil 


el>o 


"hostility" of environment 


C^Tvar 


G [0, 1] 


Geographic variability of the temperature ^ 




G [0, 1] 


global mutation rate 



Figure 2: A summary of the system parameters, which are explained in detail 
in sec. 12.31 

(b) In each individual cell, the fire lasts for one timestep and is extin- 
guished afterwards (bij —¥ 0). A burning cell which has been burning 
already in the last step is set back to uninhabited state. 

The system parameters are summarized in figure [2] 

3 Simulation and Results 

3.1 Simulation and Data Analysis 

The result of a typical simulation run is displayed in figure [3j After a simulation 
run, the following data is available: 

• The time series of carbon C in the atmosphere, 

• the time series of average temperature adaptation 

A{t) := 1 - \ nj (t) - T 3 {t)\ , (12) 

ij 

• the time series of average versatility 

n*):=j^EM*)°fc(*)' ( 13 ) 

ij 

where the defintion N^ft) :— nas been used. From these time series, 

the authors have extracted maximum, average and minimum values for A and 
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admax: 0.82244 sigmax: 10153 12 

adav: 0.7576 sigav: 0.99699 10 

admin: 0.59124 sigmin: 0.97933 g 

Pmax= 0.030383 Mpos = 9 £i 

MPDiff = 15 AvDitf = 0.15B27 4 

CJii = 0.074221 alpha = O.0GG741 2 



C.I 




50 100 150 200 250 



50 100 150 200 250 



Figure 3: Results of a simulation run for the parameter choice Ni = 12, Nj = 24, 



^max = 256, p s tart = 0.9, afi rc = 0.002, a 



growth 



= 0.1, h 



hostil 



= 1, "Tvar = 0.5, 



fi = 0.1. 

The graphical output contains the following elements: 

In the first row the system parameters are summarized, the first plot displays 
the currently stored amount of carbon, the second one shows currently inhabited 
(red/light) and uninhabited (blue/dark) cells and the third one is a graphical 
representation of the last forest fires. 

In the second row one has the Fourier spectrum (absolute values), including a 



power-law fit according to equation ( 14 ) , an illustration of the total amount of 



carbon stored in the cells over the whole simulation run and the time series of 
carbon in the atmosphere. 

The third row contains the current r-values, the current er-values and a time 



series of average temperature adaptation ( 12 ) 



The fourth row displays the time series of row-averaged r and a and a time 
series of average versatility (JT3J) . (In case that row-averaged values were not 
well-defined for certain times [when no cell of that particular row was inhabited] , 
this is reflected by white spots in the corresponding graphs, cf. figure. [8j) 



7 



V and certain characteristics of the Fourier spectrum of C (examining only 
absolute values, i.e. the amplitude spectrum X(n), which is closely related to 
the power spectrum). 

These characteristics are the position and the value of the dominant maxi- 
mum beyond n = 1 and the parameters of a power-law fit 

C 

P(n) « — + o offsot . (14) 
n a 

In addition the authors have analyzed the number of exceptional increases 
of C (usually caused by forest fires), which are defined as timesteps r, where 
C(t) > C(t — 1), but C(t — 1) < C(t — 2). and the average over the changes 
AC(r) :=C{t)-C{t-1). 

3.2 Results 

Comparison of simulation runs for various sets of parameters revealed that the 
most crucial parameters are ag re and /ihostii- Thus the authors have performed 
simulations for values of these two parameters chosen on a grid. For a wide 
mesh size and a^ vaI — \, the results are shown in figure |4j 

Some of the most important features of the model are already clearly recog- 
nizable: For afi re = there are no forest fires, so for small /ihostii all carbon is 
stored in the cells. For large values of ^hostii, it is still possible that the whole 
population dies out due to lack of adaptation. This happens in general (for the 
chosen set of parameters) for /ihostii > 1.75, regardless of afi re . 

For /ihostii = and afi rc small, but nonzero, one finds a local maximum 
of atmospheric carbon. This peak is a key feature of the model and will be 
examined closely in the following analysis. The basic explanation is that very 
small afire implies long periods without fires, in which large amounts of carbon 
are stored, such that one forest fire can have devastating consequences and 
eliminate almost all cells at once, releasing large amounts of carbon. 

The authors find adaptation maxima at /ihostii * 1 due to the fact that 
for smaller values of /ihostii the punishment for bad adaptation is less severe, 
while for even larger values of /ihostii the population more frequently dies out 
completely without having time to adapt. 

To study these effects more closely, the authors have restricted the grid to 
a smaller region in parameter space. Performing again several simulation runs 
gave the results depicted in figure [5] In addition they changed the parameter 
«Tvar to zero (figure |6| and one (figure [7]) and find that the former case (where 
temperature is solely dominated by C) is significantly different, while the latter 
case gives results very similar to those obtained for axvar = h- 
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a{C) 

q,V ■■ 



2 in: 



f(Pmax). 
■0.02, 



a(Mpos). 

.100, 



ttf" 



1 ^h^^a? 1 W 



cr(o) <j(A) 

:■ ..0.2 

L ., > * o 

at 5 - 1 ^T^-^g^ 12 ^ 




Figure 4: Results of simulation runs (10 runs for each specific combination of pa- 



rameters) for the parameter choice iV; 



12, Ni 



24, U 



256, p s . 



().' 



Qgrowth = 0.1, QTvar = 0.5, \i = 0.1. The parameters <xi = afl rc and hh = /ihostn have 
been varied as aa vc — 0.005 k, k = 0, . . . , 20 and /ihostn = 0.25 £, £ = 0, . . . 8. On the 
left hand side, the main plot shows the average carbon content C of the atmosphere. 
On the right hand side, the first row of graphs gives the amplitude of the dominant 
Fourier mode, X max , its position M poB and the exponent a of the spectral power-law 



fit (14 1. The second row shows the average temperature adaptation (A), the total 



number iVcinc of timesteps for which C increased (which corresponds to the number 
of significant forest fires) and the average change ACi nc of carbon content during such 
events. Plots of the same quantities are given one more time (second graph on the 
left hand side, third and forth row on the right hand side), but usually shown from 
a different viewing position. In the last row the standard deviations a of the seven 
quantities under consideration are given. 
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Figure 5: Results of a sequence of 10 simulation runs. Parameters are chosen 
as in figure |4j but only the most interesting region is examined by choosing 
a Rrc = 0.001 k, k — 1, ... ,20 and /ihostii = 0.1 £, t = 0, . . . 15. Again the graphs 
are displayed for two different viewing positions. 

Note that for larger values of aa IC the carbon content C is almost constant 
for small /ihostu and increases monotonically as a function of ft-hostii, starting at 
^hostii ~ 1- For small values of afi re and /ihostii, there are large fluctuations in 
C, to be noticed also in the graph for <r(C). 
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Figure 6: Results of a sequence of simulation runs. Parameters are chosen as 
in figure [5] except for axvar = 0. Again the main graphs are displayed for two 
different viewing positions. 

For this choice of parameters, the temperature does not depend on geographical 
latitude, but only on C . Apparently in this situation, for small values of ag re 
it is more likely that the atmosperical carbon level remains low, i.e. the cells 
inhabitance is stable. This is also reflected in in the Fourier amplitude spectrum, 
where the maximum still has low intensity, indicating the absence of significant 
oscillations as they can be recognized in figure [3] 
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Figure 7: Results of a sequence of simulation runs. Parameters are chosen as in 
figure [5] except for QTvar = 1. Again the graphs are displayed for two different viewing 
positions. In this case, where temperature is independent of the carbon content C, 
one finds results not so different from those displayed in figure [5] This, together 
with the significantly different results shown in figure [6] seems to indicate that the 
system characteristics change signifcantly for small values of ciTvar, i.e. the mostly 
C-dominated case. 
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3.3 Interplay of parameters 

The model parameters are not completely indepent. For example it is clear that 
in figure[4]the population typically dies out for fohostil > 1-75. One could expect, 
however, that a larger mutation rate [i may be beneficial in dealing with hostile 
environments and that one could find values of /i that allow the population to 
survive in most cases. 

This is indeed the case, as illustrated in figure [HJ For /ihostil = 1-75 and 
\i = 0.3 the population typically survives. This is also true for /i = 0.5, though 
in this case it repeatedly comes close to extinction. This is a hint that also in 
this model a mutation rate which is too large can threaten the survival of a 
species. 

4 Possible Extensions of the Model 

There are several possible extensions of this model, which might be interesting 
to study: 

• One could modify the grid by of introduction of uninhabitable spots, dis- 
tributed randomly or as clusters of some fractal shape. One may also con- 
sider the growth of such uninhabitable spots (which can be interpreted as 
spreading of cities and deserts). 

• In a more complete model one could explicitely consider the influence of 
soil, which can store carbon from deceased plants. Another important 
factor is humidity, which can reduce the rate of forest fires and also has 
influence on temperature. This could help to make the model more real- 
istic, but at the same time harder to analyze. 

• One could dismiss the normalization condition C < 1 and increase the 
total amount of Carbon available in the system during the simulation run 
(to mimic the effect of burning fossile fuels). 

• One could try to explicitely include the effects of water vapor by setting up 
a a separate cycle for H 2 0, coupled to the plant life grid, but endowed with 
a separate dynamics on a shorter timescale. It might even be worthwhile 
to think about coupling a Stochastic Cellular Automata Carbon climate 
model as proposed in this article to a more traditional model based on 
atmosphere and ocean dynamics. 
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Figure 8: Large values for /ihostii can be partially compensated by increasing the 
mutation rate /i. These plots show typical simulation runs for /ihostii = 1-75 and 
fi = 0.1, fi = 0.3 and fj, — 0.5. (Other parameters were chosen as in figure [3j) 
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• A fascinating extension may be the coupling of the system to an "ocean" , 
i.e. a carbon sink with limited absorption rate and temperature-dependent 
capacity. 

This may be especially important, since the solubility of CO2 in water 
decreases with increasing temperatur^] So while the oceans can store 
enormous amounts of CO2, their capacity decreases with global warming. 
This could lead to a dramatic positive feedback chain: 

In such a scenario, an increase of temperatures above a certain critical 
level could lead to a massive release of CO2 from the oceans and thus a 
further rise of temperature. Confronted with this frightening perspective, 
even simple models which allow to study the dynamics of such catastrophic 
behaviour can be extremely valuable. 

Apart from these extensions, there are also several aspects of the original model 
that should be studied in greater depth. The interplay of model parameters, as 
sketched in subsection |3.3[ certainly deserves closer examination. In addition it 
may be interesting to study size-dependence, in particular the correspondence 
between "critical" a n re and the extension of the grid (for example by considering 
the limit iVjJVj — > 00, a^ c — > with the product NiNja^e held fixed). 

5 Summary and Conclusions 

The authors have studied a "small-world" model designed to describe the in- 
terplay between global climate, carbon storage and sudden outbreaks of forest 
fires. 

They have seen that the least predictable behaviour is found for small, but 
nonzero values for the probability ag re of a cell catching fire. While in most cases 
an increased hostility of the environment leads to a smaller average population 
and thus to a higher carbon content of the atmosphere, this is not necessarily 
true in the most complex (and probably most realistic) case of small, but nonzero 

«firc- 

This "critical" behaviour is more dramatic in the case of temperature dom- 
inated by atmospheric carbon (instead of geographic parameters), in terms of 
this model for axvar ~ 0. 

They have also seen that there exists an interesting interplay of parameters, 
as for example the fact that a larger hostility of the environment can be par- 
tially compensated by a higer mutation rate. In general, the model exhibits a 
rich structure worth further investigation, and offers the possibility of several 
interesting extensions. 



2 This is in general true for the solubility of gases in liquids, while it is exactly the other 
way round for the solubility of solids in liquids. This is easy to understand, since enhanced 
thermal movement of the molecules makes it easier for gases to evaporate, while it makes it 
harder for solids to crystallize. 
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